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Abstract We review progress in the development of physically realistic three di- 
mensional simulated models of the galaxy. We consider the scales from star forming 
molecular clouds to the full spiral disc. Models are computed using hydrodynamic 
(HD) or magnetohydrodynamic (MHD) equations and may include cosmic ray or 
tracer particles. The range of dynamical scales between the full galaxy structure and 
the turbulent scales of supernova (SN) explosions and even cloud collapse to form 
stars, make it impossible with current computing tools and resources to resolve all 
of these in one model. We therefore consider a hierarchy of models and how they 
can be related to enhance our understanding of the complete galaxy. 

Keywords ISM • simulations • magnetohydrodynamic • MHD • hydrodynamic • 
supernova • cosmic rays • magnetic field • galaxy 

PACS 98.62.-g 

Mathematics Subject Classification (2000) 85.06 



1 Introduction 

Modelling numerically any astrophysical phenomena in general has proven chal- 
lenging, even taking into consideration some considerable simplifications. How- 
ever with the improvements in computing from the late nineties onwards reason- 
able 3-dimensional approximations of galactic features have been developed. 

The aim of this section of the chapter is to briefly review progress and to sum- 
marise some key findings to date. To impose some finite limit on the scope of the 
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review we restrict our attention to three dimensional models, which include real- 
istic parameters on galactic scales. Please attribute any omissions as a reflection of 
the ignorance of this author and not on the value of the work. 

We present a roughly chronological review of models and their methods, objec- 



tives and results. This approach evolves out of the earlier 2D study by Rosen et al. 
1 1993[ l which investigated the hydrodynamical turbulent structure of the ISM using 
stellar heating and star formation 



Tomisaka l |1998) attempted to model multiple supernovae (SNe) in 3D using a 



super-bubble approximation of in a stratified ISM with magnetic field. We sum- 
marise in Table[l] some key features of the models, which followed subsequently. 



2 Review of galaxy simulations in 3D 

In the following sections we summarise the work of several groups, who have pub- 
lished in some cases several papers over the last decade. For brevity we have cited 
the more recent publications, whilst including material from earlier papers. Please 
refer to the cited papers for the full list of citations. 



2.1 Korpi 



Korpi et al. l |1999| l constructed a comprehensive 3D model of a section of the galactic 



disc using parameters typical of the solar neighbourhood. The simulation included 
non-ideal MHD applied to a Cartesian grid of 0.5 x 0.5 x 2 kpc 3 . The resolution was 
about 8 3 parsecs 3 (pc). 

A vertical gravitational field (Kuijken and Gilmore ( 1989) ) and a shearing box 



to replicate the differential rotation were added. Thus a vertically stratified den- 
sity profile for the interstellar medium (ISM) was applied. Heating was applied via 
Supernovae (SNe) Type I and II distributed randomly in time and space. Thermal 
conductivity, kinetic and magnetic viscosities were modelled and hyper viscosity 
applied to handle shocks. 

The model replicated a multi-phase medium with the warm and hot phases in 
dynamic pressure equilibrium. The cold phase was not well reproduced through 
lack of a thermally unstable cooling function and low resolution. 

The resultant ISM structure was highly differentiated in Z. Volume filling factors 
for temperature varied with height. Correlations for mass and velocity varied with 
temperature. The outward flow of the galactic fountain appeared to be described 
within the model, but the size of the model was insufficient to recover the global 
fountain behaviour. 

In this study the warm phase was characterized by strong turbulent flow, but in 
the hot phase flows were on the scale of the remnant size enclosing the hot regions. 
This could be a feature of the numerical limits of the model. 

With this model the group also studied the evolution of super-bubbles (SBs) 
of SNe clusters and how this was affected by the strength of the magnetic field 
and turbulence. They found that the frequency of SB break outs from the disc to 
the halo was significantly increased in a turbulent regime compared to previous 
experiments in a homogeneous medium. 
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The role of magnetic field strength in confining the bubbles was reduced. This 
could be arise from the increased random element to the field due to turbulence. 
Thus our expectation of the mixing between the disc and the halo would be higher. 



2.2 de Avillez 
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Fig. 1 de Avillez and Breitschwerdt 
(black), 200 < T10 3 -" (orange), 10 s * 



2007 Scatter plot from study (2005) of B versus p for T200 
< T10 4 - 2 (blue), 10 4 ' 2 < TIO 5 - 5 K (green), and T > 10 5 5 K 



(red) regimes at 300 Myr of disc evolution. The points in the plot are sampled at intervals of four 
points in each direction. Note that during the evolution of the system the field strength broadened 
its distribution spanning two orders of magnitude after 300 Myr. 



Using a model since 2000 |de Avillez and Breitschwerdt| |2007| applied 3-D HD 
equations to a Cartesian grid 1 kpc z by 8 kpc. They used adaptive mesh refinement 
(AMR) with a resolution ranging between 1.25 3 and 10 3 pc 3 . Also included was a 
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vertical gravitational field and Type Ib+Ic and Type II SNe. 60% were clustered in 
OB associations. Thermal conductivity and kinetic viscosity were neglected, but a 
thermally stable radiative cooling function was included. 

The model evolved from hydrostatic equilibrium over a period of about 65 Myr 
to a dynamically steady state. Distinct cold (T < 10 3 K), warm (10 3 < T < 10 4 K), 
warm ionized (10 4 < T < 10 5 K) and hot (T > 10 5 K) phases were co-existing with 
varying scale heights. The cold phase was restricted to thin irregular strips around 
the mid-plane. Hot gas dominated above 1.5 kpc and warm gas in the disc below 
500 pc. 

The model was extended vertically to ±10 kpc between 2001 and 2007. They 
identified bubbles and super-bubbles typically up to 120 pc across and 200 pc high 
breaking out of the disc periodically. The cold sheet like structures that evolved, and 
which resembled the fragments of SN remnants, were in fact found not necessarily 
to be correlated with the bubble remnants, but in fact appeared to be a more general 
effect of the turbulent compressions and convection currents. 

For a further adaption of the model to MHD an initial uniform 5.8fiG magnetic 
field in the azimuthal direction was included. The horizontal magnetized disc only 
briefly delayed the disc-halo cycle observed in the HD simulations and a hot phase 
volume filling factor of between 17 to 21% resulted both with HD and MHD. 

Supersonic and super- Alfvenic flows led to strong MHD shocks. A highly tur- 
bulent field evolved with a field strength ranging between 0.1 and WpiG, which was 
not strongly correlated to the density fluctuations. This is illustrated in figure 2.2 



Magnetic pressure dominated the coldest gas, ram pressure dominant above 200 K 
and thermal pressure above 10 5 5 K. In cold regions the relation B <x p" did not hold, 
with a in the range -0.006 to 0.085. 



2.3 Balsara 



Bals araetal.| ( |2004 l applied 3D ideal MHD to a (200 pc) 3 re gion with SN randomly 
distributed. They included cooling (Raymon d et al.|[l976) and a constant diffuse 



heating term. Initial equilibrium was achieved by balancing heating and cooling for 
some given isothermal uniform density. The cooling considered could reproduce 
hot and warm phases, but not the cold phase. 

They explored the results when varying the initial temperature and density and 
applying SN rates 1 and 100 x the rate in the solar neighbourhood. They investi- 
gated resolution at about 1.5 3 and 0.75 3 pc 3 . Stratification and differential rotation 
were not considered nor thermal conductivity and kinetic viscosity. 

They found evidence of a fast dynamo sensitive to the SN rate. Over a 40 Myr 
simulation the field was amplified by two orders of magnitude. The growth rate 
was stronger in the high resolution, due to additional amplification available at 
smaller scales. 

Taking a turbulent scale from an SN remnant to be 30pc they compared energy 
growth at scales above and below this. They found that eventually the growth rates 
were independent of scale. 

Growth rates however were sensitive to SN rates, increasing up to a threshold 
where the magnetic energy was quenched. The growth of the field depended on 
the generation of helicity from the SN shocks interacting with the warm ISM. This 
was suppressed as the ISM became saturated by hot remnants. This effect had a 
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critical rate somewhere between 12 and 40 times the galactic rate. There was also 
some increased amplification if the SN rate was more intermittent and field growth 
was most vigorous where the ISM had lower temperature and higher density. 



2.4 Hanasz 
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2009). Unmagnetized regions are 



Fig. 2 Distribution of toroidal magnetic field from 1 
white, while positive and negative toroidal fields are shaded red and blue, respectively. Note that 
the color scale in magnetic field maps is saturated to enhance weaker magnetic field structures 
in the disc peripherals. The maximum magnetic field strengths are 5.9 x 10~ 4 ,4.4 x 10~ 3 , 1.5 and 
29fiG at f = 0.02, 0.7, 2.5 and 4.8 Gyr respectively. 



Form 2002 Hanasz et al. (12009) investigated dynamo action with rigid rotation 
on a magnetized galactic disc. A uniform azimuthal field was applied to a rotating 
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Cartesian domain 600 x 1800 x 600 pc with resolution 5 x 20 x 2.5 pc. They solved 
isothermal resistive MHD equations (neglecting kinetic viscosity). 

In 2004 they included cosmic rays (CR) randomly injected at spherical locations 
representing SN remnants. Differential rotation was now included in a 0.5 x 1 x 1.2 
kpc 3 Cartesian grid, with a resolution of 10 3 pc 3 . The gravitational potential was as 
described by Ferriere l |2001) and the cosmic ray diffusion tensor that of Ryu et al. 
( 2503) 

They observed an exponential growth in the magnetic field over a period of 
about 1000 Myr. Growth in the vertical component dominated, indicating that the 
CR buoyancy was important to the dynamics. Growth in the azimuthal direction 
reached 1400 x the seed field amplitude and 400 x in the radial direction. 

Growth in the ideal regime onset rapidly after 250 Myr, due to magnetorota- 
tional instability (MRI). They argued that MRI amplification is suppressed when 
comparing the resistive MHD to ideal MHD due to magnetic reconnection in the 
thin current sheets generated by compressions in the field. 

They observed that the fastest growing dynamo model resulted from periodic 
switching on and off the SN activity. Increased SN activity enhanced vertical field 
lines, increasing the loss of CR to the halo. Reduced activity supported their hori- 
zontal realignment thus enabling the amplitude of the field to grow. So the irregu- 
larity of the SN distribution can add to the dynamo. 

Hanasz et al.| ( 2509) investigated the source of the galactic magnetic seed field 



using their CR-driven dynamo model. Their domain was 50 2 x 10 kpc 3 with reso- 
lution 100 3 pc 3 . They started with an isothermal, vertically and radially stratified 
density, with a gravity potential and rotational shear, neglecting the central bulge. 

Into an unmagnetized ISM they randomly exploded SN remnants. Some 10% of 
these included a randomly oriented weak dipolar magnetic field and CR equivalent 
to 10% of the typical kinetic energy of an SN remnant. Over a 4.8 Gyr simulation 
they produced a magnetic field, which became efficiently ordered within 2 Gyr, 
aligned along the spirals and including field reversals. The vertical field evolves 
an x-shaped configuration. In their movies, vertical outflows could be seen, which 
transported the counter polarity or helicity out of the system (Figure 2.4 small in- 
serts). This seemed to be crucial for dynamo action. 

The number and orientation of the spiral configurations is sensitive to the CR 
diffusion coefficients. They could not be derived from first principles without smaller 
scales, so the assumed value for diffusivity was large and unrealistic by necessity. 
Snap shots of the configuration are shown in Figure 2.4 



2.5 Slyz 



In their 2002 2D model Slyz et al. l |2005| l investigated the characteristics of the spiral 
arms through isothermal HD. They used a 23 2 kpc 2 with a resolution of 115 2 pc 2 , 
including kinematic viscosity and applying a centrifugal and gravitational poten- 
tial. 

They varied the temperature, via the sound speed, the density profiles and rates 
of galactic rotation. Given the typical sound speed throughout the galaxy is far 
exceeded by the angular velocity of the gas relative to the spiral arms, it might 
not have been expected to have much effect the structure. They found that the 
non-axisymmetric pattern of the spiral galaxies diminished with increasing sound 
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speed. At levels above 25 to 30 km s they faded away. As sound speeds reduced 
the mass in the spiral arms increased. 

They also found that the inner galaxy was more sensitive to changes in the frac- 
tion of galactic mass in the disc. As the disc fraction increased, velocities at the 
centre increased massively whilst the structure in the outer regions was fairly in- 
sensitive to the disc fraction. 

Subsequently Slyz et al. (2005) conducted a more localized analysis of star for- 
mation and SN feedback with a 3D hydrodynamical model. They used a 1.28 3 kpc 3 
Cartesian grid together with a particle mesh to track the stellar mass. The model 
included a thermally unstable cooling function and self-gravity. Stratification and 
any external gravity were not included. 

A random Gaussian velocity perturbation was applied at the start onto an ISM 
of uniform density and temperature. They found that the regime saturated to a 
multi-phase medium, with velocity dispersion generally settling to levels consistent 
with the temperature of the medium: 15 km s _1 , 30 and 75 for cold (< 2000 K), 
warm (< 10 5 K) and hot (< 4 x 10 6 K) respectively. Velocities in the very hot (> 
4 x 10 6 K) are far more dispersed, with up to 500 km s _1 observed. These violent 
flows strip away elements of the cold dense star forming regions as they interact. 

The distribution of mass was also a feature of the ISM temperature. The bulk 
of the volume was filled with hot diffuse medium of below 10~ 3 atoms cm 3 . The 
bulk of the mass was restricted to dense cold filamentary or sheet-like or worm-like 
structures accounting for a very small portion of the volume. In the absence of SN 
shocks the probability density functions (PDF) for the ISM were well described by 
the log normal. However with SN feedback a much higher incidence of increased 
density regions were observed, suggesting that star formation rates might be en- 
hanced in a turbulent ISM. 



2.6 Mac Low 



Continuing the work from Balsara et al. l |2004[ l, Mac Low et al. l |2005} studied the 



pressure distributions in the turbulent ISM with slightly differing results to those 
of Slyz et al. 12005). They applied ideal MHD equations in a (200 pc) 3 domain, with 
3.12, 1.56 and 0.83 parsec resolution, including a radiative cooling function and uni- 
form uv-heating, but excluding gravity. 

They randomly injected SN remnants, at 1 and 4 times the solar neighbour- 
hood rate and analysed the evolution of the pressure distribution. They found the 
distribution of the pressure anomalies was independent of the orientation of the 
magnetic field due to the super-Alfvenic flows. Typically the hot gas formed dis- 
crete clumps, which were enclosed by warm or cold gas. The distribution of the 
colder dense clusters, tended not to be closely aligned directly to remnant struc- 
tures. These formed on scales of dozens of parsecs and appeared to result from the 
general effects of large scale turbulence. 

The presence of dense isolated filaments was more effective at breaking the mo- 
mentum of remnants than equivalent mass more smoothly distributed and even 
with 4 times the SN rate hot zones remained unconnected and warm gas remained 
the largest volume filling factor. Even with SN driving they found the PDFs of the 
pressure were well described by the log normal and not a power law. The disper- 
sion range was an order of magnitude away from the mean, far from equilibrium. 
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2.7 Piontek 

The turbulence in the galaxy is not adequately accounted for by SN activity alone so 
in 2005 Piontek and Ostriker ( 2007) modelled MRI-driven turbulence in the warm/ cold 
phases of the ISM. They apply ideal MHD to a (200pc) 3 domain with resolution of 
0.78 and 1.56 pc and realistic heating and cooling. They model differential rotation 
for the solar neighbourhood. The growth of turbulence was studied in a medium 
in uniform density and pressure and also one in which randomly distributed cold 
spherical clouds were included amongst the warm ISM. The MRI instability rapidly 
dominated and the saturated state was independent of the initial state. They stud- 
ied the relationship between the MRI driven turbulence and the two-phase thermal 
structure derived from the cooling function by comparison to MPJ is a single phase 
environment. 

They found that the gas separated into two phases with peaks at 100 K and 
8000 K, with an increasing spread of gas lying outside the equilibrium as turbulence 
saturated. The relative proportions of gas in the two phases was dependent on mean 
density, with the proportion of mass in the cold phase increasing substantially as the 
mean density increased. 

Pressures were in reasonable equilibrium even in two phases without turbu- 
lence, but in the saturated state the pressure dispersion increased five fold. 

Magnetic field growth 0.26fiG to about 2.5^G was not sensitive to density nor 
phase. Away from SN activity MRI could certainly be a factor in the turbulent 
structure. Extrapolating their results for mean densities typical in the inner galaxy, 
R < 10, kpc they would predict turbulent velocity scales of about 4kms~ 1 , which 
could be a significant component. 

Piontek and Ostriker l |2007} included a vertical gravity profile, which now in- 



duced a vertical stratification in both density and temperature. The turbulent mix- 
ing of the warm ISM had a significantly higher scale height than the cold, which 
now was more clustered than in the earlier model. This combined with MRI to in- 
crease the magnetic pressure and they argued that the MRI could be responsible 
for suppressing star formation in the less dense outer region of the mid plane, even 
given the presence of cold clouds. 



2.8 Joung 



From 2006 Ryan J oung et al.| ( 2009) conducted a HD study in a vertically stratified 



ISM including stellar and halo gravity and realistic diffuse heating and cooling plus 
SN turbulence. The vertical profile of the heating was realistically motivated. They 
neglected thermal conduction and kinetic viscosity, self-gravity and differential ro- 
tation. 

They found the scale height of the galactic fountain to be several kpc. The vol- 
ume filling factor near the mid-plane of the hot gas was inversely related to the 
magnitude of the diffuse heating, which they tested over a range up to 4 times that 
of | Wolfire efaT] fl995) . This was due to a shift in the thermal equilibrium position 



to exclude the cold phase and increase the volume component of the warm gas. 

The vertical density profile was out of agreement with observations which they 
attributed to lack of pressure, in the absence of magnetic and cosmic ray pressure. 
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In this study and also from 2009 they find that neither the density nor the kinetic 
energy wavelengths follow the Kolmogorov power law, but instead are far more 
widely dispersed. 



2.9 Gressel 
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Fig. 3 Vertical slices from abov e the mid-plane (u pper panels) and horizontal mid-plane slices (lower 
panels) after t = 112 Myr from\Gressel et al.\i2008\. From left to right [in logarithmic scales]: number 
density[-4.82,1.10]cm~ 3 , column density I17.34,21.55]cm~ 2 , temperature[2.06,7.20]K, velocity dispersionl- 
0.15,2.77]km s _1 & magnetic field strength[-4.16,0.34]uG. 



Gressel et al. ( 2008 1 conducted a non-ideal 3D MHD simulation using the NIR- 
VANTA code. Their domain was 0.8 2 x 4 kpc 3 and included realistic cooling, thermal 
conductivity and kinematic viscosity and rotational shear. A snapshot from their 
simulation is shown if figure 2.9 to illustrate the resolution and structure obtained 
through such simulations. 

They tested a variety of SN rates (1/4 to 1 x uq) and rotational frequency (1 to 
8 x £lo), where ctq and Oo denote the observed rates in the solar neighbourhood, 
and measured the resulting dynamo parameters. 

A result common to each run, was that the velocity dispersion increased with 
height up to about 1 kpc, even though the peak of SN activity was much lower. They 
found that turbulent pumping was directed inward, but balanced by a outward 
wind of comparable magnitude. 
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They found the dynamo not to be sensitive to the SN rates considered, 1, 2 & 4x 
the observed rates. However the dynamo increased only with galaxy rotation rates 
2 & 4x the observed rate. They did not rule out however that these conclusions 
could also be sensitive to the assumed gravity and density profiles. They found no 
dynamo in the models without rotational shear. 



3 Summary 



At the highest level o f resolution we have had the study of Balsara et al. |2004|, 



Piontek and Ostriker 



(2007) and 



Mac Low et al. 



1 2005 ' in 200 J pc 3 domains. Too 



small to include large scale structure, such as density stratification or spiral arms, 
these studies were able to relate the relationship of turbulence from SN or MRI, 
with dynamo amplification. 

They were able conclude that MRI could be a significant source of turbulence 
in the galaxy and in particular in the outer galaxy away from SN activity. The pres- 
ence of cold dense filaments, they demonstrated, were not closely correlated to the 
dynamics of the SN remnants. Instead they appeared to be a product of generalised 
turbulence and gravitational instability. 

Finally they found dynamo amplification rate was sensitive to the level SN rate 
up to an upper bound of between 12 & 40 x the solar neighbourhood rate. Ampli- 
fication of the magnetic field was also increased by greater intermittency in the SN 
rate. Independent of SN rate the dynamo is quenched once the ISM is saturated by 
hot 

JT999),|de Avillez and Bre- 



Losing resolution, but gaining structure Korpi et al, 
|itschwerdt| l |2007| , |Slyz et al.| ( |20()5) , |Ryan Joung et a 



(20091 



and Gressel et al. 



( |2008) constructed models incorporating the density stratification of the ISM near 
the plane of the galactic disc. 

Unlike the smaller models, the ISM does not become saturated by SN activity. 
The stratified structure permits the disc to expand and relax as the SN rate fluctu- 
ates in time and spatial distribution. This relieves the quenching effect of turbulent 
saturation on the magnetic dynamo. 

However results for the dynamo have been mixed. Korpi et aL"| (1999 1 and Gres- 
sel et al.| {2008} did not find a dynamo with solar neighbourhood model parame- 



ters, and |Gressel et al.| ( |2008) only found the dynamo with differential rotation 2 
and 4 x the galactic rate Both these models used non-ideal MHD, where magnetic 
reconnection may dampen the dynamo in the critical cold dynamic filaments. The 
numerical values for resistivity are of course much higher than the realistic values 
observed in the ISM. de Avillez and Breitschwerdt| l |2007} did not investigate the 
dynamo but included a 5/ muG magnetic field using ideal MHD within their later 
work. 

With a vertical range of 20 kpclde Avillez and Breitschwerdt|l|2007| were able to 
include elements of the galactic fountain and without boundary mass losses, could 
sustain simulations over 400 Myr. They found that an excess of 200 Myr was re- 
quired to completely saturate and establish an equilibrium state for the disc-halo 
cycle. They identified the scale height of the disc-halo interface at about ±1.5 kpc. 
The additional height was required to allow the hot gas in the halo to cool and rain 
back to the disc. 
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Despite outflow boundary conditions, and hence mass losses over the duration 
of the simulation, Gres sel et al.|p008) were able to sustain their simulations over 1 
Gyr. Once the turbulent state had evolved they found the boundary flows were min- 
imal, so much of the dynamical cycle could be contained within a domain height of 
±2 kpc. 

On the grand scale Hanasz et al. l |2009) were able to produce a familiar model 
of the spiral structure of the galaxy. Restricted to isothermal models, they found the 
structure of the spirals sensitive to the CR diffusion. In addition the 2D study of 
Slyz et al. (2005) found the structure was also dependent on the temperature of the 



ISM, as parameterized by sound speed. 

These models have informed our understanding of how the ISM behaves on 
different scales across a range of parameters. It is also clear that the interaction of 
temperature, density, magnetic field and cosmic rays combined are significantly 
altered in the absence of any one of these with unexpected consequences. 

The challenge is to successfully incorporate all these elements over a number 
of physically meaningful scales. Global models of galaxies are indispensable, with 
sufficiently high resolution to include the small-scale physics. To resolve SNe rem- 
nants requires a resolution under 5 3 pc 3 . For the random field dynamo we may re- 
quire scales < 2 3 pc 3 . The gravitational field (and self -gravitation of the gas) needs 
to be included, the latter requiring scales below lpc. Density waves are expected 
to interact with the magnetic fields and need to be considered. Arm and inter-arm 
modelling would require domains extending a few kpc in the galactic plane. To in- 
clude even the lowest reach of the halo, would require a vertical range of at least 4 
kpc. 
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